As already demonstrated by experiment,^[44] KIIIIK peptides capped with an acetyl in the N terminal and an amine in the C terminal can self-assemble to form monolayer nanotubes, so the intra- and inter-sheet structures shown in Fig. 1 have been constructed as the initial configurations for studying the intermolecular forces involved in the nanostructures by means of MD simulations. The initial intra-sheet structure composed of two peptides was constructed as perfectly aligned parallel β-sheet structure, which is the most probable conformation appeared in the replica exchange MD simulation of a KIIIIK trimer.^[45] The initial inter-sheet structure was constructed by face-to-face stacking of two cross-β sheets, which contains six strands perfectly aligned to form parallel β-sheets. All these initial structures were constructed using the PyMol^[46] and VMD^[47] programs.

All the simulations were performed with the GROMACS simulation package.^[48] The KIIIIK peptides with protonated lysine were modelled by the OPLS-AA force field,^[49,50] and the water molecules were modelled by the TIP4P force field.^[51] A rectangular box with periodic boundary conditions in all three dimensions was applied and a 2-fs time step was used for all the simulations. The cutoff for the short-range non-bonded interactions was 1.2 nm, and the long-range electrostatic interactions were treated by the particle mesh Ewald algorithm.^[52,53] Chloride ions modeled by the OPLS-AA force field were added to the solution to neutralize the system. The NPT ensemble was applied with its temperature maintained by a Nose-Hoover thermostat^[54,55] and its pressure maintained by a Parrinello–Rahman barostat.^[56] To check the stability of the two structures, 2-ns simulations with unconstrained peptide positions were conducted after a 0.1-ns simulation with constrained peptide positions to allow the water molecules to relax around the peptides.

The last configuration of the 2-ns simulation served as the initial structure of the subsequent steered MD simulations, in which the peptides should be placed in such a position that there is an adequate space for the pulling simulation along the reaction coordinate. As shown in Fig. 1, in the steered MD simulations, both intra- and inter-sheet structures had a fixed part and a moving part, with the latter pulled with a constant rate. The pulling direction for the intra-sheet structure is along the H-bond direction, which resides on the vector from the center-of-mass (CoM) of the fixed strand to the CoM of the pulled strand, as illustrated in Fig. 1(a). The pulling direction for the inter-sheet structure is along the side-chain direction, which resides on the vector from the CoM of the fixed sheet to the CoM of the pulled sheet, as illustrated in Fig. 1(b). The pulling rate of 5 nm/ns and the spring constant of 1000 kJ/(mol/nm²) were applied to the 0.5-ns steered MD simulations for calculating PMFs. To study the mechanical response, simulations with a pulling rate of 0.25 nm/ns were also conducted for each structure. The final CoM distances between the fixed and pulled parts were approximately 3.0 nm for the intra-sheet structure and 3.5 nm for the inter-sheet structure.

On the basis of the non-equilibrium steered MD simulations, the umbrella sampling technique^[39,40] and the weighted histogram analysis method (WHAM)^[57] were used to calculate the free energy profiles (i.e., PMFs) along the reaction coordinates. The initial configuration in each window of the umbrella sampling was selected from the steered MD trajectory. Totally 26 windows were used in the umbrella sampling calculation for the intra-sheet structure with the CoM distance ranging from 0.44 nm to 3.24 nm, and each window has a width of approximately 0.108 nm to allow adjacent windows to have an adequate overlap. There were totally 23 windows for the inter-sheet structure ranging from 1.27 nm to 3.15 nm, and the window width was approximately 0.082 nm. Inside each window, a 5-ns NPT simulation with the simulation parameters described above was performed.

To obtain the intermolecular bond stiffness, the lowest potential wells of the PMFs around the equilibrated structures were fitted by a harmonic function, 1 where k is the elastic constant corresponding to the intermolecular bond stiffness, ξ is the reaction coordinate defined as the CoM distance between the pulled part and the fixed part, ξ₀ is the equilibrated CoM distance between those two parts, and G₀ is the free energy minimum located at the equilibrium conformation.

Since the NPT ensemble was adopted in the simulations, the PMFs actually depict the profile of the one-dimensional Gibbs free energy G along the reaction coordinate. A Δ symbol will be added before each thermodynamic variable in the following equations because the free energy obtained from umbrella sampling is a relative value. The entropy along the reaction coordinate can be calculated by the finite difference of the temperature derivative of the free energy^[58] 2 where ΔG is the Gibbs free energy, ξ is the reaction coordinate, T is the temperature, and ΔT is the temperature difference between two simulations. In this work, both the intra- and inter-sheet structures were simulated at two temperatures: 285 K and 300 K, respectively. The enthalpy contribution can then be obtained by 3

The loading rate of the external force significantly affects the peak force of the force profile in a steered MD, which is the minimal force required to break the structure.^[35,37] To minimize the non-equilibrium effect, we conducted another set of steered MD simulations with a very slow pulling rate of 0.25 nm/ns to study the mechanical response. The results shown in Figure 2 demonstrates three stages for the dissociations of both structures under an external pull at a constant rate. For the intra-sheet structure, from 0 ns to 1 ns the force linearly increases up to a peak value but ξ has little increase; from 1 ns to 3.5 ns the force decreases in a complex way and ξ increases linearly with a rate of 3.6 Å/ns; from 3.5 ns to 10 ns the force fluctuates around zero and ξ increases linearly with a rate of 3.6 Å/ns. For the inter-sheet structure, from 0 ns to 2 ns the force increases linearly to a peak value but ξ has little increase; from 2 ns to 5.2 ns the force decreases in a complex way and ξ increases linearly with a rate of 3.5 Å/ns; from 5.2 ns to 10 ns the force fluctuates around zero and ξ increases linearly with a rate of 3.8 Å/ns. As shown in Fig. 2(a), the peak force for the intra-sheet structure is about 211 kJ/(mol/nm), a little smaller than about 300 kJ/(mol/nm) for the inter-sheet structure, and the times required to reach the peak value are about 1 ns and 2 ns, respectively. Since the peak force equals the minimum force required to break the structure, and the intra-sheet structure contains only one pair of peptides but the inter-sheet structure contains six pairs of peptides, the minimal force required to break an individual pair of the intra-sheet structure along the H-bonding direction is 211 kJ/(mol/nm) and the corresponding minimal force for the inter-sheet structure along the side-chain direction is 50 kJ/(mol/nm), indicating that the intermolecular forces in the hierarchical self-assembled structures are anisotropic, qualitatively agree with the experimental results obtained by Fitzpatrick et al.^[27]

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Pulling forces and CoM distances as a function of time in the two steered MD simulations with a pulling rate of 0.25 nm/ns at T = 300 K.

The PMFs for the intra- and inter-sheet structures have been calculated and are shown in Fig. 3. For both the intra- and inter-sheet structures, there exist almost no free energy barriers during peptide association along the H-bond and side-chain direction. The free energy differences between the associated and dissociated states are −65.3 kJ/mol and −62.4 kJ/mol for the intra- and inter-sheet structures, respectively.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Total PMFs at T = 300 K and harmonic fittings of the potential wells for the intra- (a) and inter-sheet (b) structures.

Each PMF in Fig. 3 and S1 in Appendix A exhibits a potential well marked by two dashed lines, which are located at around 0.46 nm for the intra-sheet structure and 1.2 nm for the inter-sheet structure. As shown by the red lines in Fig. 3 and S1, the potential wells can be well fitted by using Eq. (1), so we can estimate the optimal conformation parameters (the CoM distance and the free energy for the equilibrated structures) and the intermolecular bond stiffness for both structures at 285 K and 300 K. All the fitted optimal conformation parameters and intermolecular bond stiffness at 300 K are listed in Table 1, along with the experimental results for the intra- and inter-sheet structures of an amyloid peptide VQIVYK,^[27] since no experimental results are available for KIIIIK. The optimal conformation parameters ξ₀ for both structures obtained from the simulations agree well with experimental values. The simulation values of the intermolecular bond stiffness k are apparently larger than the experimental values, but both have the same high ratio of about twenty for the intra-sheet value over the inter-sheet value, demonstrating that the intermolecular forces for the two structures are highly anisotropic, consistent with both experimental and simulation results.^[27,35] The quantitative difference between our simulation results and the experimental results may be attributed to the difference in peptide types studied by simulation and experiment as well as the fact that peptides in simulation are solvated in water whereas peptides in experiment are dry samples.

Table 1.

Table 1.

Table 1. Intermolecular bond stiffness values, optimal conformation parameters, and minimum free energies. The data in brackets are experimental results for peptide VQIVYK.[27] . Parameters Inter-sheet Intra-sheet k/(N/m) 0.32 (0.07) 5.58 (2.0) ξ0/nm 1.27 (0.80–1.20) 0.45 (0.47) G0/(kJ/mol) −65.3 −62.4

Table 1. Intermolecular bond stiffness values, optimal conformation parameters, and minimum free energies. The data in brackets are experimental results for peptide VQIVYK.[27] .

According to the thermodynamics relation in Eq. (2), the entropy contribution ΔΔS₁₂ to free energy difference between two states can be calculated by For a spontaneous transition from state ξ₁ to state ξ₂, if ΔΔS₁₂ > 0, the entropy is the driving force for the association; if ΔΔ S₁₂ < 0, the enthalpy is the driving force. To understand the thermodynamic features of the driving forces for the intra- and inter-sheet structures, we have conducted umbrella sampling simulations for the intra- and inter-sheet structures at the temperatures of 285 K and 300 K, and the PMFs for the two structures are shown in Figs. 4(a) and 4(b). According to the PMFs for the intra-sheet structure in Fig. 4(a), the free energy difference between the associated and dissociated states is −49.3 kJ/mol at 285 K, larger than −65.3 kJ/mol at 300 K. The intermolecular bond stiffness value of the intra-sheet structure at 285 K obtained by fitting the potential well is 6.08 N/m, which is larger than the value of 5.58 N/m at 300 K. The decrease of the free energy difference between the associated and dissociated states with increasing temperature demonstrates that entropy difference between the two states is positive, indicating that entropy contributes to the driving force for peptide association along the H-bond direction. The decrease of the intermolecular bond stiffness with increasing temperature implies that the free energy difference between the equilibrated structure and the state around equilibrated structure increases with increasing temperature, indicating that enthalpy is the driving force keeping the intra-sheet structure stable. The PMFs for the inter-sheet structure in Fig. 4(b) indicate that the free energy difference between the associated and dissociated states is −54.6 kJ/mol at 285 K, which is also larger than −62.4 kJ/mol at 300 K. The intermolecular bond stiffness value of the inter-sheet structure at 285 K is 0.09 N/m, which is smaller than the value of 0.32 N/m at 300 K. The decrease of the free energy difference between the associated and dissociated states with increasing temperature demonstrates that entropy difference between the two states is positive, indicating that entropy contributes to the driving force for peptide association along the side-chain direction. The increase of the intermolecular bond stiffness with increasing temperature implies that the free energy difference between the equilibrated structure and the state around equilibrated structure decreases with increasing temperature, indicating that entropy also contributes to the driving force keeping the inter-sheet structure stable.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) PMFs at T = 285 K and 300 K for the intra-sheet structure (a) and the inter-sheet structure (b), as well as the entropy and enthalpy contributions at T = 300 K for the intra-sheet structure (c) and the inter-sheet structure (d).

To identify the thermodynamic driving forces for the association processes along the reaction coordinates in detail, the entropy (−TΔS) and enthalpy (ΔH) contributions for both structures in the association process were obtained by Eqs. (2) and (3) based on the PMFs at 285 K and 300 K. Since the PMFs for the associated states are smaller than those for the dissociated states, when −TΔS decreases with decreasing ξ, the association is driven by entropy; when ΔH decreases with decreasing ξ, the association is driven by enthalpy. From Figs. 4(c) and 4(d), we can see that the curves of the entropy and enthalpy contributions for both structures undulate, indicating that enthalpy and entropy in turn drive the association of the strands for the intra-sheet structure along the H-bond direction and the cross-β sheets for the inter-sheet structure along the side chain direction. The results for the formation of intra-sheet structure in Fig. 4(c) show that −TΔS decreases but ΔH increases when ξ decreases from 0.82 nm to 0.7 nm, indicating that the association from 0.82 nm to 0.7 nm is entropy driven. From 0.70 nm to 0.54 nm, −TΔS increases but ΔH decreases, indicating that the driving force is enthalpy. From 0.54 nm to 0.47 nm, −TΔS decreases but ΔH increases, indicating that the driving force is entropy. From 0.47 nm to 0.38 nm, −TΔS increases but ΔH decreases, indicating that the driving force is enthalpy. Comparing −TΔS and ΔH at the associated state when ξ is about 0.45 nm with those at the dissociated state where ξ is larger than 2.5 nm, we can see that −TΔS for the associated state is smaller than that for the dissociated state, whereas ΔH is larger. These results indicate that entropy overall drives the whole association process, but enthalpy stabilizes the equilibrium intra-sheet structure at the distance of about 0.45 nm. Similar to the discussion for the formation of intra-sheet structure, the association of cross-β sheets can be divided into five ranges whose driving forces are entropy, enthalpy, entropy, enthalpy, and entropy in the ranges from 2.13 nm to 1.97 nm, 1.97 nm to 1.78 nm, 1.78 nm to 1.70 nm, 1.70 nm to 1.52 nm, and 1.52 nm to 1.29 nm, as shown in Fig. 4(d). Comparing −TΔS and ΔH at the associated state for ξ around 1.27 nm with those at the dissociated state with ξ larger than 2.5 nm, we can see that −TΔS for the associated state is smaller than that for the dissociated state, but ΔH is larger. These results indicate that entropy also drives the formation of the inter-sheet structure for the whole association process, and entropy stabilizes the equilibrium inter-sheet structure at about 1.27 nm.

The VDW and electrostatic interactions between different parts have been analyzed in detail to unravel the anisotropy of the intermolecular bond stiffness and the difference in thermodynamic driving forces for the intra- and inter-sheet structures. The short-range pairwise interactions between peptide–peptide E_PP, peptide–solvent E_PS, and solvent–solvent E_SS along the association reaction coordinates have been calculated and are shown in Figs. A2–A4 in Appendix A. All the energy differences between the associated and dissociated states are clearly different, as listed in Table 2. Here the short-range non-bonded interactions contain the VDW and electrostatic energies between atom pairs whose distances are less than 1.2 nm. The E_a energy of the associated state is defined as the average energy of structures with ξ smaller than 0.5 nm for the intra-sheet structure or smaller than 1.4 nm for the inter-sheet structure. The E_d energy of the dissociated state is defined as the average energy of structures with ξ larger than 2.5 nm for both structures. From Table 2 we can see that, for both structures, ΔE_PP (peptide–peptide) and ΔE_SS (solvent–solvent) are negative, but ΔE_PS (peptide–solvent) is positive, indicating that the total short-range interactions between peptides promote association and only the peptide-solvent interaction promotes dissociation.

Table 2.

Table 2.

Table 2. Short-range non-bonded energies of the associated and dissociated states and their differences. . Intra-sheet Inter-sheet (kJ/mol) Ea Ed ΔE Ea Ed ΔE P−P VDW −134 0 −134 −191 0 −191 Coul −121 0 −121 −104 0 −104 P−S VDW −252 −377 125 −766 −1002 236 Coul −1354 −1640 286 −7293 −7497 204 S−S VDW 83884 83865 19 243106 243076 30 Coul −526573 −526372 −201 −1524433 −1523991 −442 Total −444550 −443924 −36 −1288711 −1289414 −281

Table 2. Short-range non-bonded energies of the associated and dissociated states and their differences. .

To determine the contributions of the hydrophobic side chains and the hydrophilic backbones of KIIIIK to E_PP, the interactions between side chains, E_ScSc, and those between backbones, E_BbBb, are now calculated. Starting with the intra-sheet structure, we can see that, with increasing ξ increases towards zero whereas decreases towards zero, as shown in Fig. 5(a). The energy differences between the associated and dissociated states and are −25 kJ/mol and 2.5 kJ/mol, respectively, indicating that interactions between hydrophobic side chains promote the formation of the intra-sheet structure through the VDW interaction. For the interaction between hydrophilic backbones, both and increase towards zero as ξ increases (Fig. 5(c)) and and are −30 kJ/mol and −115 kJ/mol, respectively, indicating that the interactions between hydrophilic backbones are dominated by the electrostatic interactions and promote the formation of the intra-sheet structure. Since is significantly smaller than , the electrostatic interactions between hydrophilic backbones should dominate the formation of the intra-sheet structure.

Fig. 5.

Figure Option ViewDownloadNew Window

Fig. 5. (color online) Short-range non-bonded energies between different groups of the pulled and fixed parts as a function of ξ for the intra- and inter-sheet structures at T = 300 K. Panels (a) and (c) are non-bonded energies of the hydrophobic side chain-side chain and hydrophilic backbone-backbone interactions for the intra-sheet structure, and panels (b) and (d) are those for the inter-sheet structure.

For the inter-sheet structure, Figure 5(b) shows that, with increasing ξ, increases towards zero but changes little from around zero. is about −110 kJ/mol, indicating that the interactions between hydrophobic side chains promote the formation of the inter-sheet structure and mainly come from the VDW interactions. For the interactions between hydrophilic backbones, both and show little change from around zero as ξ increases, indicating that the interactions between hydrophilic backbones have little effect on the formation of the inter-sheet structure.

The competition between the peptide-peptide and peptide-solvent interactions has a significant influence on the peptide conformations and association of peptides in the hierarchical processes.^[59,60] The ensemble average of the short-range non-bonded energy of the solvent in the shell around the peptide from x nm to x + 0.1 nm will be denoted as E_PSx. According to E_PSx from 0.2 nm to 0.6 nm shown in Figs. A5–A8 in Appendix A, we find that both the VDW and electrostatic energies between peptide and solvent in different shells are sensitive to ξ and variations of these non-bonded energies with respect to ξ are different in different shells. To distinguish the roles of the hydrophilic backbones and the hydrophobic side chains in the peptide-solvent interactions, the short-range non-bonded energies between hydrophobic side chains and solvent in shells from 0.2 nm to 0.6 nm (E_ScSx) and those between hydrophilic backbones and solvent in the shells (E_BbSx) were calculated for each umbrella-sampling window.

For the hydrophobic side chains in the formation of the intra-sheet structure, increases but , , and decrease with ξ, as shown in Fig. 6. The energy differences between associated and dissociated states are −11, 15, 10, and 13 kJ/mol for , , , and , respectively. These results indicate that only the VDW interaction of the hydrophobic side chain with water in the shell from 0.2 nm to 0.3 nm promotes the formation of the intra-sheet structure, but those with water in the other shells promote the dissociation of the intra-sheet structure. The electrostatic energies, , , , and show little change with ξ in Fig. 6, indicating that the electrostatic interaction of hydrophobic side chains with water has little effect on the formation of the intra-sheet structure.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Short-range non-bonded energies between the hydrophobic side chains and solvents around the strands in different shells as a function of ξ for the intra-sheet structure at T = 300 K. Panels (a), (b), (c), and (d) correspond to shells from 0.2 nm to 0.3 nm, 0.3 nm to 0.4 nm, 0.4 nm to 0.5 nm, and 0.5 nm to 0.6 nm, respectively.

For the hydrophilic backbone in the formation of the intra-sheet structure, increases but , , and decrease with ξ, as shown in Fig. 7. The energy difference between the associated and dissociated states are −38, 5, 25, and 11 kJ/mol for , , , and , respectively. These results indicate that the VDW interaction of the hydrophilic backbones with water in the shell from 0.2 nm to 0.3 nm promotes the formation of the intra-sheet structure, but those with water in the other shells promote the dissociation of the intra-sheet structure. The electrostatic energies in Fig. 7, , , , and decrease with ξ and the electrostatic energy differences between associated and dissociated states are 20, 46, 65, and 58 kJ/mol for , , , and , respectively. These results indicate that the electrostatic interactions of hydrophilic backbones with water in all the shells promote the dissociation of the intra-sheet structure.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) Short-range non-bonded energies between the hydrophilic backbones and solvents around the strands in different shells as a function of ξ for the intra-sheet structure at T = 300 K. Panels (a), (b), (c), and (d) correspond to shells from 0.2 nm to 0.3 nm, 0.3 nm to 0.4 nm, 0.4 nm to 0.5 nm, and 0.5 nm to 0.6 nm, respectively.

For the hydrophobic side chain in the formation of the inter-sheet structure, increases but , , and decrease with ξ as shown in Fig. 8. The VDW energy differences between associated and dissociated states are −17, 60, 63, and 36 kJ/mol for , , , and , respectively. These results indicate that only the VDW interaction of the hydrophobic side chain with water in the shell from 0.2 nm to 0.3 nm promotes the formation of the inter-sheet structure, but those with water in the other shells promote the dissociation of the inter-sheet structure. For the electrostatic energies, , , , and fluctuate around zero, as shown in Fig. 8, and show little change as ξ increases, indicating that the electrostatic interactions of hydrophobic side chains with water have little effect on the formation of the inter-sheet structure.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. (color online) Short-range non-bonded energies between the hydrophobic side chains and solvents around the sheets in different shells as a function of ξ for the inter-sheet structure at T = 300 K. Panels (a), (b), (c), and (d) correspond to the shells from 0.2 nm to 0.3 nm, 0.3 nm to 0.4 nm, 0.4 nm to 0.5 nm, and 0.5 nm to 0.6 nm, respectively.

For the hydrophilic backbone-solvent in the inter-sheet structure in Fig. 9, all the interaction energies show little change with ξ and fluctuate around 234, −68, −163, −102, −64, −288, −394, and −322 kJ/mol for , , , , , , , and , respectively. These results indicate that both the VDW and electrostatic interactions of the hydrophilic backbones with water have little effect on the formation of the inter-sheet structure.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. (color online) Short-range non-bonded energies between the hydrophilic backbones and solvents around the sheets in different shells as a function of ξ for the inter-sheet structure at T = 300 K. Panels (a), (b), (c), and (d) correspond to shells from 0.2 nm to 0.3 nm, 0.3 nm to 0.4 nm, 0.4 nm to 0.5 nm, and 0.5 nm to 0.6 nm, respectively.